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We calculate the momentum diffusion coefficient for heavy quarks in SU(3) gluon plasma at 
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extracted from a Monte Carlo calculation of the correlation function of color electric fields, in the 
leading order of expansion in heavy quark mass. Systematics of the calculation are examined, 
and compared with perturbtion theory and other estimates. 
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1. Introduction 

The charm and bottom quarks play a very important role as probes of the medium created in 
relativistic heavy ion collision experiments. Since the masses of both these quarks are considerably 
larger than the temperatures estimated to have been reached in RHIC and even in LHC, one expects 
these quarks to be produced largely in the early state of the collision. They therefore allow us to 
look at the medium at its early times. Besides, since the nature of energy loss and of collective 
behavior of the heavy quarks are different from those of the light quarks, a study of the heavy 
quark jets and the flow of the heavy hadrons leads to crucial insights into the way the medium 
interacts. 

As collision with a thermal quark does not change the energy of a heavy quark substantially, 
one would expect the thermalization time of the heavy quark to be large. Since elliptic flow devel- 
ops early, it is reasonable to expect that the elliptic flow will show a mass hierarchy: v\ 3> 3> vf , 
where D and B refer to generic mesons in D (one valence charm) and B (one valence bottom) 
family, and h refers to the light hadrons. Similarly, the nuclear suppression factor can be intuitively 
expected to be closer to 1 for the heavy-light mesons: R 1 ^ <C R^ <C R^. 

Since the typical energy loss in a hard collision with the thermal particles is ~ T, for a thermal 
heavy quark with M S> T, p ~ VMT, it takes a large number of collisions for the heavy quark 
to change its momentum by ^(1). Therefore, one can picture scattering with thermal quarks as 
uncorrected momentum kicks, and use a Langevin description for the motion of the heavy quark 
in the thermal medium []l|, ||] : 

^ = -voPi + m, mm')) = ^ s( t - t '). q.d 

The momentum diffusion coefficient, fc, is related to the correlation of the force term: 

1 f°° 

K = 3 ]_J* L<W(0)>. (1-2) 

fc has been determined in perturbation theory^. The fluctuation-dissipation theorem relates r]o 
andK[|]: 

y\ D = -^—. (1.3) 
2MT 

The heavy quark elliptic flow, and Raa, has been measured in RHIC [Q]. While a Langevin 



formalism based description, Eq. (|1.1|), seems to describe the pj dependence of the heavy quark 
flow quite well, the value of K needed to describe the data [|J] is much larger than the leading 
order (LO) perturbation theory prediction. Recently, the next-to-leading order (NLO) calculation 
of K has been performed in the static quark limit j^]. While the NLO result is, encouragingly, 
much larger, the large change between LO and NLO also points to the unreliability of perturbation 
theory in the temperature regime of a few times the transition temperature T c . A non-perturbative 
evaluation, if possible, will greatly add to our understanding of the medium response to a heavy 
quark probe. 

Here we report the results of an non-perturbative estimation of K using lattice QCD in the 
quenched approximation (i.e., for a gluon plasma). In the next section we outline the methodol- 



2 



Heavy quark momentum diffusion coefficient 



Saumen Datta 



ogy, and in the following section we discuss the results. A more detailed discussion, including 
examination of various systematic errors, can be found in Ref. 

2. Calculational Details: 



The calculation of the heavy quark momentum diffusion coefficient, fc, is a nontrivial problem 
for lattice QCD. On the lattice, one only calculates the Euclidean-time Matsubara correlators, while 
K, and other transport coefficients, are directly connected to the real-time retarded correlators of 
suitable currents [f7|]. An analytical continuation of the Matsubara correlator to real time is required 
to extract K from it. Such a continuation of numerical data is very difficult (see [Q] for a recent 
review). On top of that, K is related to the width of a narrow peak in the spectral function of the 
quark number current, making it very difficult to extract it reliably. 

A strategy of extracting fc for an infinitely heavy ("static") quark has been formulated in Ref. 
[B, 0] which alleviates this second problem somewhat. In the static limit, the propagation of heavy 



quarks is replaced by Wilson lines, and the correlator of Eq. ( |1.2| ) reduces to the evaluation of re- 
tarded correlator of electric fields connected by Wilson lines [Jsj] . This correlator can be analytically 
continued to Euclidean time ft The lattice discretization of the Euclidean correlator leads to 
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tf(P,T)E i (T ) 0)tf(T,0)E/(0,0) 



(2.1) 



where U(t\, T2), the timelike gauge connection between points (3c, Ti) and (3c, T2), is the phase factor 
associated with the evolution of an infinitely heavy quark in imaginary time. Ej(z,x) is the color 
electric field at point (3c, t) and L = tr £/ (/3 , 0) is the Polyakov loop. 

In order to calculate G\ at {x), we use the discretization of the electric field ^ 



E t (3c, t) = Ui (3c, t) Ua (3c + i, t) - t/ 4 (3c, t) £/; (3c + 4) 



(2.2) 



which has good ultraviolet properties. We calculated G^ at (T) on a set of SU(3) pure gauge lattices 
in the temperature range 1 — 2T C . It is imperative to use sufficiently fine lattice spacings, a, so that 
we get a large number of points, N t , in the % direction. We could reliably extract K only from 
lattices with N t = 20 or more. The different lattices for which we will be quoting estimates of fc 
are shown in Table [jj 
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N, 


T/T c 


# sublattices 


# sublattice updates 


6.76 


20 


1.04 


5 


4000 


6.80 


20 


1.09 


5 


3000 


6.90 


20 


1.24 


5 


2000 


7.192 


24 


1.50 


4 


2000 


7.255 


20 


1.96 


5 


2000 



Table 1: List of lattices which were used to calculate K from Ge(x). Also given are the parameters used 
for multilevel update: the number of sublattices the T direction was divided in, and the number of sublattice 
averagings before a full lattice update. 
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A precise calculation of G^ at (r) on lattices with such large N t is known to be very difficult 
with a naive updating algorithm. We adapted the multilevel algorithm which was devised 



precisely for such problems, to calculate Eq. < \2.l\ ). The lattice is divided into several sublattices 



The expectation value of the correlation functions are first calculated in each sublattice by averag- 
ing over a large number of sweeps in that sublattice while keeping the boundary fixed. A single 
measurement is obtained by multiplying the intermediate expectation values appropriately. The 
number of sublattices and the number of sublattice averagings were tuned for the various sets, so as 
to get correlators with a few per cent level accuracy. The parameters used for the lattices in Table |T] 
are also shown in the table. The use of the multilevel algorithm turned out to be absolutely essen- 
tial for our calculation: for large z, it was up to ^(10 3 ) times more efficient than a naive updating 
algorithm. For details about implementation of the algorithm, and its performance, see Ref. [B]. 
Before extracting K, we need to convert G\ at {z) to the physical correlator of the electric field, 

G E {z)=Z{a)GY{z) (2.3) 

where Z(d) = Z| is the lattice spacing dependent renormalization factor for the electric field cor- 



relator. We use the tree-level tadpole factor for Z(a) Ql 1Q. A nonperturbative evaluation of Z(a) is 
planned for the future. 

3. Results: 

Standard fluctuation-dissipation relations ^ connect the momentum diffusion coefficient, K, 
to the low-ft) behavior of the spectral function: 

dco , . coshwfT- 1/2T) 

7" ( " sl# (31) 

2T 

K = lim — pica). (3.2) 
co->0 CO 

The nontrivial part of the calculation is to extract p(ffl) from Ge(z). In this work we assumed 
a functional form for p(ft)), so that calculation of p(co) and K become a fitting problem. We use 
the simple fit form 

p(ffl) = |ffl0(A-ffl) + bco\ (3.3) 



where the first term in the r.h.s. of Eq. ( J3-3| ) is the low-ft) diffusion part, and A is an infrared cutoff. 
As elaborated later, we fix A = 3T for our central results and fit for K, b. For the fit, x 2 minimization 
was carried out with the full covariance matrix included in the definition of % 2 . We typically 
obtained acceptable fits to the correlators for %a in the range [N t /4,N t /2], with # 2 /d.o.f ~ 1. At 
shorter distances, lattice artifacts start contributing and the simple form of Eq. ( |3.3[ ) does not work 
well. Also using the leading order lattice correlator instead of the continuum form did not improve 
the quality of the fit. We, therefore, restrict ourselves to the long distance part of the correlator. 

In order to get a feel for the contribution of the diffusive part of the spectral function to the 
correlator, in Fig. |l] a) we show the correlators reconstructed from different parts of p(ft)) sepa- 
rately. In this figure we take the N t = 24, l.5T c data set, and use the best fit form of Eq. ( |3.3[ ) (for 
A = 3T). The contributions to the total correlation function from the ft) 3 part of p(ft)) and that 
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Figure 1: (a) Ge{t) and the correlator reconstructed from the diffusive part(DIFF), shown normalized to 
that reconstructed from the leading order part (LOC), at 1 .5 T c (j3 = 7. 192); see text, (b) The contribution of 
the diffusive part, normalized by the leading order part, at the different temperatures in Table [|. 



from the diffusive part are calculated separately using Eq. (3.1). In the figure these two parts are 
referred to as LOC and DIFF, respectively, and the total correlator and DIFF are shown normalized 
by LOC. While the correlator is dominated by the contribution from the boo 3, term, the diffusion 
term has a substantial contribution near the center of the lattice. In Fig. [T] a) it contributes ~ 17 % 
at xT = 0.5. 

In Fig. [I] b) we show the ratio of the diffusive and the leading order parts of the correlator with 
varying % at the different temperatures studied in Table [jj Here we show the 1-a bands and not the 
best fit values. At all temperatures, except the one at the highest temperature, the diffusive part is 
seen to reach about 5 % level by xT ~ 0.3. Note that the accuracy of our correlator is better than 
this. Also no significant trend of temperature dependence is seen in this figure, indicating the lack 
of a strong temperature dependence of K in this temperature range. 

The value of the momentum diffusion coefficient, k, obtained from the analysis outlined 
above, is shown in Fig. || a). The central error bar corresponds to the purely statistical error, 
obtained from a Jackknife analysis. The larger error bar corresponds to uncertainties due to various 
systematics: 

• As mentioned above, for the central value shown in the figure, we have set A = 3T. In 
order to look for the dependence of the result on this choice, we have varied A in the range 
[2!T,oo). A systematic error band which envelops the central values of the fits with varying A 
is introduced. 



The fit form, Eq. (p.3|), is a simple model, taking into account the leading order perturbative 
behavior and the fact that at small ft), p((o) °< ft). In order to test the dependence of the 



extracted value of K on Eq. ( |33| ) we also use an alternate form for the infrared part of p (ft)): 

p 2 (ft>) = fctanh^ &(co-A) + boo 3 . (3.4) 
This form is motivated by classical lattice gauge theory [|12|]. The difference between K 



obtained from uses of Eq. (p.3|) and Eq. (3.4) is also included in the systematic error. 
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Figure 2: (a) k/T^ extracted from Gjb(t), as explained in the text, at various temperatures in the gluon 



plasma, (b) The diffusion coefficient, D, obtained from K using Eq. (3.5 ). Also shown are the LO perturba 



tion theory estimate [|J| and a different lattice estimate jl5[|, both for gluon plasma. The experimental value 
quoted by PHENIX, and an estimate from the jV = 4 supersymmetric Yang-Mills theory are also shown. 



• For the central values in Fie. 0, we have used the range T m j n to where T m j n is the 

smallest % for which we got a good % 2 . We looked for the stability of the fit values for 
variation of T m j n . In all sets except one, we could get stable fits, with good % 2 for T m j n > N r /4. 
For the set at 1.96 T c , where the fit value was not so stable, we included the variation of the 
fit value into our systematic error estimation. 

A detailed discussion of the sizes of the various systematics can be found in Ref. Within 
the small variation of LT our resources permitted, we did not find a finite volume dependence above 
our other errors for LT > 2. 

The value of k/T at 1.5 T c , shown in Fig. ^|a), agrees within errors with a similar calculation 
by Francis et al. [13], while an earlier calculation Jl4| ] found smaller values. It is also an order-of- 
magnitude larger than the LO perturbation theory value that can be extracted from Ref. [0L 

The experimental results for heavy quark diffusion are usually presented in terms of the dif- 
fusion coefficient D, which controls diffusion of the heavy quark in position space. The Einstein 
relation, 

T IT 2 

D = —— = , (3.5) 

Mr\ D k 

connects D with k. D obtained from Eq. ( |3~5| ) and Fig. |2| a) is shown in Fig. || b). The leading 
order perturbation theory value |Q] is also shown there. 

A more direct approach to calculating D from lattice QCD would be to look at the correlation 
function of the heavy quark number current, QjiQ. A calculation of D for the gluon plasma using 
such an approach has been presented in Ref. Jl5|]. Their results are also shown in Fig. || b), 
where we have added there statistical and systematic errors in quadrature. The results are even 
further from the perturbation theory estimate and systematically lower than our results, although 
consistent within the large error bars. 

Our results are for a gluon plasma, so a comparison with the experimental results [Q] requires 
due caution. At the least, a comparison of the results in Fig. ^| with the perturbative results for 
quenched QCD gives us an indication of how much the nonperturbative results can change from 
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the perturbative results in the deconfined plasma at moderate temperatures < 2T C . Even then, the 
results are most encouraging since they indicate that the nonperturbative estimate for DT can easily 
be an order of magnitude lower than perturbation theory, bringing it in the right ballpark to explain 
the V2 data. 

Since dimensionless ratios of various quantities are known to scale nicely between quenched 
and full QCD if plotted as function of T /T c , one can, more optimistically, hope that our results, 
as plotted in Fig. ||, will be even quantitatively close to the full QCD values. In this spirit, in Fig. 
^| b) we also compare the lattice results with the experimental data. The lattice results seem to be 
a little above the best range for description of the PHENIX data [Q] using the Langevin approach 
[0], though reasonably close within our large systematics. Interestingly, our lattice results seem to 
show very little temperature dependence in the temperature regime studied here. 

As mentioned in Sec. [IJ, the heavy quark diffusion coefficient has also been calculated in a 
very different theory, the jV = 4 supersymmetric Yang-Mills theory with the number of colors, 
N c — > °°, at large 't Hooft coupling = g 2 N c , using AdS/CFT correspondence [@J. Of course, 
this theory is very different from QCD in many respects. Moreover, it crucially exploits symmetries 
which QCD does not have. However, to get a feel for what kind of values of D such a theory would 
predict for parameters relevant for QCD at ~ l.5T c , we naively set N c = 3 and as = 0.23 in the 
results of Ref. [|p. This gives DT ~ 0.2 (Fig. ^|b), which is lower than, but in the same ballpark as 
our estimate. 
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